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A criterion for the choice of optimal softening length e and time-step dt for A-body 
simulations of a collisionless stellar system is analyzed. Plummer and Hernquist spheres 
are used as models to follow how changes in various parameters of an initially equilibrium 
stable model depend on e and dt. These dependences are used to derive a criterion for 
choosing e and dt. The resulting criterion is compared to Merritt's criterion for choosing 
the softening length, which is based on minimizing the mean irregular force acting on a 
particle with unit mass. Our criterion for choosing e and dt indicate that e must be a 
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! factor of 1.5-2 smaller than the mean distance between particles in the densest regions to 
be resolved. The time-step must always be adjusted to the chosen e (the particle must, on 
O ■ average, travel a distance smaller than 0.5e during one time-step). An algorithm for solving 
N-body problems with adaptive variations of the softening length is discussed in connection 
with the task of choosing e, but is found not to be promising. 
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1 Introduction 

The evolution of gravitating systems on time scale shorter than two-body relaxation time 
is described by the collisionless Boltzmann equation. The number of particles in iV-body 
simulations of such systems is usually several orders of magnitude smaller than the number of 
stars in real systems. It could be said that in such N-body simulations collisionless Boltzman 
equation is solved by Monte-Carlo method, in which the particles are used to represent mass 
distribution of real system. In practice, integrating the equations of motion of the particles 
involves smoothing the acutal (pointwise) potential of each individual particle, for example, 
by substituting a Plummer sphere potential for the actual potential. This procedure modifies 
the law of gravity at small distances: 

F *f = Gmmj Jl^i- _> Flf = Gm imj ru r i~ ri , 2W2 , (1) 
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where F£J and F-° — are the real and softened forces, respectively, acting on a particle of 
mass rrii, located at point r i; produced by a particle of mass rrij, located at point Tj, and e 
— is the softening length determing the degree of potential smoothing. 

Potential smoothing is used in iV-body simulations for two reasons. 

First, attempts to solve the equations of motion of gravitating points using purely New- 
tonian potentials and simple constant-step integration always lead to problems during close 
pair encounters (the integration scheme diverges). Correct modeling of such systems requires 
the use of either variable-step integrating algorithms or sophisticated regularization meth- 
ods, which lead to unreasonably long CPU times. No such problems arise when a softened 
potential is used. 

Second, smoothing the potential reduces the "graininess" of the particle distribution, 
thereby making the potential of the model system more similar to that of a system with a 
smooth density distribution, i.e., a system described by the collisionless Boltzmann equation. 

It is obvious that e cannot be too large: this would result in substantial bias of the 
potential, and would also impose strong constraints on the spatial resolution of structural 
features of the system. It is important to have an objective criterion for choosing the softening 
length in iV-body simulations. In this paper, we derive such a criterion by analyzing the 
time variations of the distribution functions for spherically symmetric models with different 
values of e . 

2 Merritt's criterion for choosing e 

Merritt || proposed a criterion for choosing the softening length based on the minimization 
of the mean irregular force acting on a particle. He introduced the mean integrated square 
error (MISE), which characterizes the difference between the force produced by a discrete 
set of N bodies and the force acting in a system with a continuous density distribution: 



Here, E denotes averaging over many realizations of the system, F(r, e) is the force acting 
on a particle with unit mass located at the point r produced by N particles with softened 
potentials, and F trU e(r) is the force acting on the same particle in a system with a continuous 
density distribution p(r). 
Two factors contribute to the MISE. 

- Fluctuations of the discrete density distribution. These are very important at small 
e, decrease with e (see the left-hand branch of the curve in Fig. Q), and for a given e, 
decrease with increasing N [Hj; 

- Errors in the representation of the potential (the difference between the softened po- 
tential and a point-mass potential). These errors, on the contrary, are very large for 
large e, decrease with decreasing e, and do not depend on N (the right-hand branch of 
the curve in Fig. ^) 

The MISE reaches its minimum at some e = e m (Fig. [TJ). Merritt suggested to adopt this e m 
as the optimal choice of the softening length in TV-body simulations. This quantity depends 
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Figure 1: Dependence of the MISE on the softening length e for a Plummer model with 
N = 1000. 



on N and, e.g., in the case of a Plummer model, can be approximated by the dependence 
in the interval N from 30 to 300 000 fin virial units, where G — 1, the total 
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mass of the model is M tot = 1, and the total energy of the system is E x 



tot 



-i/4) m 



3 Formulation of the problem 

Merritt's criterion is based on minimizing the mean error in the representation of the force 
in an equilibrium system at the initial time. It is of interest to derive a criterion for choosing 
the softening length based directly on the dynamics of the system. If a stable equilibrium 
configuration is described by a discrete set of N bodies, due to various errors, the parameters 
of the system (e.g., the effective radius) will deviate from their initial values with time. The 
smaller these deviations, the more adequately the model reproduces the dynamics of the 
system. It is reasonable to suppose that these variations will be minimized for some e. We 
analyzed two equilibrium models: a Plummer sphere model 
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as an example of a model with a nearly homogeneous density distribution, and a Hernquist- 
sphere model with an isotropic velocity distribution [Hj 
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as an example of a model with more concentrated mass distribution. In (3) and (4), ap and 
an are the scale lengths of the density distributions for the Plummer and Hernquist models. 
(The Plummer and Hernquist models contain about 35% and 25% of their total masses inside 
the radii a P and an, respectively). We use here virial units (G — 1, M tot = 1, E tot = —1/4), 
for which ap = 37r/16 ~ 0.59 and an = 1/3. We computed the gravitation force using 
the TREE method In the computations, we set the parameter 8, which is responsible 
for the accuracy of the computation of the force in the adopted algorithm, equal to 0.75. 
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Test simulations made with higher (6 = 0.5) and lower (9 = 1.0) accuracy showed that 
our conclusions are independent of 9 within the indicated interval. Hernquist et al. [Til] 
showed that the errors in the representation of the interparticle force in the TREE method 
are smaller than the errors due to discreteness noise, independent of N, so that the above 
approach can be used to solve the formulated problem. We integrated the equations of 
motion of the particles using a standard leap-frog scheme with second order accuracy in 
time-step. The number of particles was iV = 10 000, and the computations were carried out 
in the NEMO [Zj software package. 

We investigated two methods for smoothing the potential. The first is based on substitut- 
ing the Plummer potential for the point-mass potential, and modifies the force of interaction 
between the particles according to the scheme (JTJ). The second method uses a cubic spline 
|Hj to smooth the potential. We report here results only for the first method, but all our 
conclusions remain unchanged in the case of spline smoothing. 

We analyzed how various global parameters of the system are conserved, in particular, 
the variation of the density distribution and, separately, the variation of the particle- velocity 
distribution. 

We characterized the deviation of the particle distribution in space at time t from the 
initial distribution using the quantity A r , which was computed follows. We subdivided the 
model into spherical layers and computed for each layer the difference between the number 
of particles at time t and the number at the initial time. We then computed A r as the sum of 
the absolute values of these differences and normalized it to N, the total number of particles 
in the system. The thickness of the layers was 0.1 and the maximum radius was equal to 2 
(each of the spherical layers in both models contained a statistically significant number of 
particles, always greater than 50). There is an important remark: A r will not be equal to 
zero for two random realizations of the same model. This quantity has appreciable natural 
noise due to the finite number of particles considered. We estimated the level of this natural 
noise to be the mean A r averaged over a large number of pairs of random realizations of the 
model. 

We computed the parameter A„ for the distribution of particles in velocity space in 
a similar way. For the results reported here, the thickness of the layer was 0.1 and the 
maximum velocity was 1.5 (this choice of maximum velocity ensured that each spherical 
layer contained a statistically significant number of particles for both models more than 30). 

We also computed the two-body relaxation time for the models as a function of e . 
We determined this times scale as the time over which the particles deviate significantly 
from their initial radial orbits (the deviation criterion and method used to estimate it were 
similar to those employed by Athanassoula et al. (3]). The estimation method can be briefly 
described as follows. We constructed a random realization of the system (the Plummer 
or Hernquist sphere). In this system, we chose a particle moving toward the center in an 
almost radial orbit and located about one effective radius from the center. Then the orbit 
of this particle was touched up to strictly radial. We followed the evolution of the entire 
system until this particle traveled a distance equal to 1.5 times its initial distance from the 
center. Knowing the time t p that has elapsed since the start of the motion and the angle of 
deflection of the particle from its initial radial orbit, a p , we can compute the time required 
for the particle to deviate from its initial direction by one radian. We then average the 
resulting time over a large number of similar particles to estimate the two-body relaxation 



4 



time 

1 _ / _Op_ 
^relax \ y/^p 

where t re iax is the two-body relaxation time and < . . . > denotes averaging over a large 
number of test particles. Formula © is valid in the diffusion approximation, i.e., assuming 
small angles for individual scattering acts. 

4 Results of numerical simulations 
4.1 The Plummer Sphere 

Figures El - 0] illustrate the results of the computations for the Plummer sphere. Our main 
conclusions are the following. 

1. Deviations of the softened potential from the Newtonian potential become important 
when e > e m . The system adjusts to the new potential by changing its distribution 
function. The adjustment is almost immediate, and occurs on time scales t « 20 
(one time unit corresponds approximately to the crossing time for the core, ap). Such 
a strongly smoothed system that has evolved far from its initial state shows hardly any 
further changes (Fig. EJ). This result agrees well with the behavior of the fractional 
error of the total energy of the system for large e, discussed by Athanassoula et al. [2] 
(Fig. 7 in their paper). 

2. Close pair encounters are computed incorrectly when the particle travels a distance 
greater than e during one time-step dt. This introduces an error, which rapidly accu- 
mulates and results in substantial changes in A r and A^ (FigEJ). This happens when 
e < edt = <r v dt, where a v is the mean velocity dispersion in the system (a v « 0.7 in the 
adopted virial units for the Plummer model). 

3. The evolution of the system outside the interval < e < e m is simulated incorrectly. 

4. A comparison of the relaxation time of the system (Fig. 0} with the time scale for the 
variations of the density distribution in the system (A r ) shows that, in the interval 
eat < e < e m , A r varies only insignificantly on time intervals a factor of two to three 
shorter than the relaxation time. When e = e m , the system preserves its structure 
for the longest time, since the relaxation time is maximum for this e (upper panel in 
Fig. ED. 

5. At the same time, it is clear from Fig. El that, for the model considered, e m exceeds 
the mean distance between particles in the central regions, so that we expect strong 
modification of the potential in a substantial fraction of the system at e = e m . As 
they adjust to the changed potential, the particles rapidly change their dynamical 
characteristics (A v ; lower panel in Fig. EI)- It follows that e m is by no means the 
best choice from the viewpoint of conserving the initial velocity distribution function, 
although A„ remains approximately constant after the adjustment has ended. 
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Figure 2: Dependence of A r and A^ on e for different times for the Plummer-sphere model. 
The horizontal line shows the natural noise levels for A r and A„. The vertical lines indicate 
the position of e m — the optimal softening length according to the criterion of Merritt (for 
N = 10 000, e m = 0.05 in virial units), and the minimum e = eat to which the time-step 
was adjusted (e dt = 0.006). Averaging was performed over several models for each e. The 
time-step is dt = 0.01. 
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Figure 3: Average distance to the nearest neighbor at various radii for the Plummer and 
Hernquist models with N = 10 000. The vertical lines correspond to the radii of spheres 
containing 0.01, 0.1, and 0.5 of the total mass of the system. The horizontal lines shows the 
optimal e = e m according to the criterion of Merritt. 



6. The effect of the two factors affecting the MISE (the "graininess" of the potential and 
its modification as a result of smoothing) are separated in time. The first factor acts 
on time scales comparable to the relaxation time (for a given e), whereas systematic 
errors in the representation of the potential appear almost immediately. It follows that 
if for the chosen e systematic errors are small (i.e. e is small enough) then the iV-body 
model very accurately represents the initial system, but only over a time interval that 
is a factor of 2-3 shorter than the relaxation time for the this e (of course its right only 
if time-step was small enough, i.e. e > 

7. Choosing e to be a factor of 1.5-2 smaller than the mean distance between particles in 
the densest region yields the optimal solution from the viewpoint of conserving both A r 
and A„. For the Plummer sphere with N = 10 000, this corresponds to e ~ 0.01 — 0.02 
(in virial units). We can see from Fig. Elthat the simulation of the system is almost 
optimal when e ~ 0.01 — 0.02, but this is achieved by restricting the evolution time of 
the system ( t < 100). When e = e m ~ 0.05, the system is simulated satisfactorily on 
time scales t < 300, but we must pay for this increased evolution time: the system is 
not correctly represented by the model from the very start. 

4.2 The Hernquist Sphere 

Figures 3-5 illustrate the computation results for the Hernquist-sphere model. Our main 
conclusions are the following. 
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Figure 4: Relaxation times for the Plummer and Herquist models as functions of e. The 
vertical line indicates the optimum e m according to the criterion of Merritt. The time-step 
was dt = 0.001. Averaging was performed over 2500 and 1000 test particles for each e for 
the Plummer and Hernquist models, respectively. 
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Figure 5: Same as Fig. [2] for the Hernquist-sphere. The time-step is dt = 0.001. In this 
case, we have for N = 10 000 and this time-step e m = 0.0087 e^t = 0.0008. 
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1. We can see from Fig. El that, with the Hernquist-sphere model, e m is comparable to 
the mean distance between particles in the very central regions of the model (which 
contain less than 1% of all particles). This means that the MISE can be affected even 
by a small number of particles with strongly modified potentials. 

2. Figure El shows that the change in the particle distributions in the configuration and 
velocity space (A r and A„) are indistinguishable from the noise even for e > e m ~ 
0.0087, right to e w 0.02. This is explained by the fact that, when e = 0.02, only 10% 
of the particles have e values greater than the mean distance between particles, and 
only this small number of particles have strongly modified potentials. We simply do 
not observe the adjustment of these central regions to the new potential. 

3. We showed in the previous section that, for nearly uniform models, e should be chosen 
about a factor of 1.5 - 2 smaller than the mean distance between particles in the 
densest regions. This criterion cannot be directly applied to significantly nonuniform 
models, since the central density peaks would be washed out due to the finite number of 
particles. For such models, e should be about a factor of 1.5 - 2 smaller than the mean 
distance between particles in the regions to be resolved. In the case of a Hernquist 
sphere with N = 10 000, the iV-body model with e = 0.01 adequately represents the 
system everywhere except for the most central regions, which contain 2-3% of all the 
particles (Fig. Ej). When e = 0.02, already the regions containing 10% of all particles 
are not simulated correctly. 

4. Unlike the Plummer sphere, in the case of the Hernquist sphere, significant changes 
in A r and A^ occur on times of the order of two "relaxation times" t = 100 (Figs. 4 
and 5). We put the term relaxation time in quotes here, because this concept is not 
fully applicable in the Hernquist model. The relaxation time is, by definition, a local 
parameter, and has diferent values in the central region and periphery for the strongly 
nonuniform systems we analyze here. 

5 Do N-body simulations need adaptive code? 

A number of authors have discussed the possibility of developing a code for iV-body simula- 
tions with an adaptively adjustable softening length; i.e., with e varying as a function of the 
mean distance between particles at a given point of the system [S], j^j. Dehnen [S] showed, 
based on the formalism of Merritt [6\, that the use of such an adjustable softening length 
decreases significantly the role of the irregular forces averaged over the entire system. It 
follows from the results reported here that the only possible advantage of such a code might 
be a reduction of computing time when modeling very inhomogeneous models. This follows 
from the fact that one can choose larger softening lengths, and thereby longer time-steps, in 
less dense regions. 

We encountered a problem when implementing such a code, however. The change of e for 
individual particle results in an asymmetric change of the total energy of the system. When 
e is increased (e sma n — > e mean ), particles located near a radius of about 2e mean pass from a 
"harder" to a "softer" potential, i.e., the absolute value of the potential energy decreases. A 
decrease in e (e^g — > e mean ) results in the opposite effect, but, in this case, on average, the 
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change in the potential influences the greater number of particles located near the radius 
2ebi g , so that the absolute value of the total potential energy of the system gradually increases 
with time. Our iV-body simulations showed significant increases in the total energy. Various 
artificial methods can be suggested to compensate for this energy change, but it remains 
unclear how this will affect the evolution of the system as a whole. 

6 Conclusions 

The results of our study indicate that the softening length e in iV-body simulations should 
be a factor of 1.5 — 2 smaller than the mean distance between particles in the densest regions 
to be resolved. In this case, the time-step must be adjusted to the chosen e (on average, the 
particle must travel a distance smaller than one-half e during one time-step). 

The use of code with variable softening lengths appears to be of limited value, and does 
not show particular promise. 
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